The dynamics of an ion chain in a harmonic potential 
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Cold ions in anisotropic harmonic potentials can form ion chains of various sizes. Here, the 
density of ions is not uniform, thus the eigenmodes are not phononic-like waves. We study chains 
of N 1 ions and evaluate analytically the long wavelength modes and the density of states 
in the short wavelength limit. These results reproduce with good approximation the dynamics 
of chains consisting of dozens of ions. Moreover, they allow to determine the critical transverse 
frequency required for the stability of the linear structure, which is found in agreement with results 
obtained by different theoretical methods [D.H.E. Dubin, Phys. Rev. Lett. 71, 2753 (1993)] and by 
numerical simulations [J. P. Schiffer, Phys. Rev. Lett. TO, 818 (1993)]. We introduce and explore 
the thermodynamic limit for the ion chain. The thermodynamic functions are found to exhibit 
deviations from extensivity. 



I. INTRODUCTION 

Coulomb crystals are organized structures of charged 
particles, which interact through the Coulomb repulsion 
and organize in regular patterns at sufficiently low tem- 
peratures in presence of a confining potential 1]. These 
potentials are realized by means of Paul or of Penning 
traps 0, and their geometry determines the crystal's 
structure. Several remarkable experiments have reported 
crystallization of ion gases in Paul trapsjj, |j EJ 
Penning traps [8| and ion-storage rings The Bragg- 
scattering in three-dimensional structures was studied 
providing information about the internal structure of the 
crystal [ToL ITU . These crystals represent a kind of rar- 
efied condensed matter, the inter-particle distance being 
of the order of several micrometers, allowing to study 
the structure by means of optical radiation. Variation 
of the potential permits to control the crystal shape as 
well as the number of ions, thus offering the unique op- 
portunity of studying the transition from few particles 
to mesoscopic systems. Besides, these structures have 
been object of growing interests as t hey provide promis- 
ing applications for spectroscopy [l3 . ll3j] . frequency stan- 
dards |l4j , study and control of chemical reactions , 
and quantum information processors 0, 0, 0, . 

In this work we investigate the dynamics of Coulomb 
chains. These are one-dimensional structures, obtained 
by means of strong transverse confinement and that usu- 
ally consist of dozens of ions localized along the trap 
axis |3, |6j . They represent a peculiar crystallized struc- 
ture: In fact, due to the axial potential the equilibrium 
charge distribution is not uniform |2Ct l2l| . This is in 
contrast to the three-dimensional case, where the den- 
sity of charges in a harmonic potential is uniform and, 
therefore, where the eigenmodes are phononic-like waves. 
In the Coulomb chain the non-uniformity of the density 
of ions combined with the long-range interaction result 
in excitations that are fundamentally different from the 
phonons in solids and lead to interesting thermodynamic 
properties. The exploration of these excitations and of 



the chain thermodynamics is the subject of the present 
paper. 

Our starting point is the ions equilibrium configura- 
tion evaluated in J2jj. We investigate the dynamics for 
small oscillations, when the harmonic approximation is 
valid, in the limit of large number of ions. We derive 
analytically the eigenfrequencies and the corresponding 
eigenmodes for the long wavelength excitations. These 
are compared with numerical results and good agreement 
is found. From the resulting dispersion relation the den- 
sity of states of the long wavelength eigenmodes is deter- 
mined. An analytic form of the density of states is also 
found for the short wavelength excitations. This result al- 
lows to evaluate the critical aspect ratio between the fre- 
quencies of the transverse and the axial confining poten- 
tial, that determines the stability of the chain. The value 
we find agrees with numerical results j^], which have 
been experimentally verified for chains of few ions |23| . 
In particular, it is in agreement with the analytical esti- 
mate in |24j , which was obtained under different require- 
ments. 

Using these results we discuss the statistical mechanics 
of the chain, and derive some thermodynamic quantities 
in a specific thermodynamic limit, that is defined here by 
keeping constant the density of ions in the chain center, 
as the number of ions tends to infinity and the axial fre- 
quency to zero, analogously to the definition for cold neu- 
tral gases in traps |25J. Non-extensive thermodynamic 
properties are found. We compare the thermodynamic 
functions with the ones of a chain of finite number of 
ions that are obtained numerically, and find reasonable 
agreement. 

This work is organized as follows. In section II we in- 
troduce the basic equations and discuss the fundamental 
properties. In section III the spectrum of excitations is 
studied. In section IV we investigate the statistical me- 
chanics of the system. Section V presents the conclusions 
and the outlook. In the appendices, several details of the 
calculations of section III are reported. 
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II. A STRING OF CHARGES IN A HARMONIC 
POTENTIAL 

The Hamiltonian describing the dynamics of a chain of 
TV ions of mass m and charge Q, which are confined by a 
harmonic potential, is given by 



N 2 



3=1 



(1) 



where rj = (xj, yj, Zj), Pj are the positions and conjugate 
momenta (j — 1, . . . , N). The term V accounts for the 
oscillator's potential and the Coulomb repulsion, 
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V&i - Xj) 2 + (yi - y 3 ) 2 + (zi - Zj) 2 

where the harmonic oscillator has rotational symmetry 
around the x-axis with axial and transverse frequencies 
v, v tl respectively. 

For sufficiently low kinetic energy crystallization oc- 
curs. The temperature, at which the gas is crystal- 
lized, corresponds to large plasma parameters T = 
Q 2 j a-wsksT 1. Here, aws is the Wigner-Seitz radius, 
that is a function of the mean density n and is defined 
as aws = (3/47T71) 1 / 3 p|. In this regime the ions are 
localized at the classical equilibrium positions ij , that 

satisfy the equations dV/dr$\ (oj = 0, and such that the 

j 

potential energy is minimal. When the harmonic poten- 
tial is sufficiently asymmetric, i.e. for v -C vt, the ions 
equilibrium positions are confined to the trap axis |22| . 



namely rj ^ = (x^ UJ ,0, 0), and satisfy the equation de- 
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where the numbering convention is Xi > Xj for i > j. 
The stability of these points with respect to the trans- 
verse vibrations depends on the number of ions N and 
on the ratio vt/v : and it is discussed in Section fill Bl 
In this section, we assume that the configuration is sta- 
ble and approximate the potential by its second order 
Taylor expansion around the points rj ^. We denote by 

Qi = Xi — xf^ the displacements in the cc-direction, and 
approximate the Hamiltonian Q as 

H 1=3 Vo + -Hhar, 

where Vq = V(rj , . . . ,rjy ) is the classical minimum 
energy, while -ffhar describes the (classical) harmonic os- 
cillations around the equilibrium points |26t 1271 128| , 
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and the coefficients Kij — d 2 V/dx 2 \^ x oy are positive and 
take the form 
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Equation (0J shows that the axial motion is decoupled 
from the transverse motion in the harmonic expansion. 
The corresponding equations of motion are 
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and describe a system of coupled oscillators, with 
long range interaction and position-dependent coupling 
strength. In this paper eigenmodes will be calculated. 
For this we assume q^t) = J e lwt g i (w)dw/27r. To simplify 
notations we replace q~i (oj) by qi . This results in equations 
for the eigenmodes of frequency u> that are similar to JHJ , 
but with q\ replaced by —uj 2 qi. The same replacement 
will be performed for yt(t) and Zi{t). 

It can be easily verified that the center-of-mass motion 
is an eigenmode of the secular equations (|6I8[) . The axial 
center-of-mass mode is q\ = . . . = qisr at the character- 
istic frequency is v, while the transverse center-of-mass 
modes are yi = ■ ■ ■ = Vn and z\ = . . . = zn at frequency 
Uf We remark that the axial and the transverse coupling 
terms appearing in Eqs. (KilSH have opposite signs. Due 
to this property, in the axial direction the collective exci- 
tations are of higher frequencies than the center-of-mass 
frequency is, while in the transverse plane the collective 
excitations are of lower frequencies than v t . 



A. Properties and symmetries 

Hamiltonian Q is not translationally invariant, and 
this is a consequence of the non-uniformity of the ions 
equilibrium distribution, due to the harmonic force ap- 
pearing in Eq. (J3J. The Hamiltonian @ is however in- 
variant under reflection with respect to the center of the 
trap, which coincides with the origin of the axes. In par- 
ticular, 

T (o) _ __(0) 

where % — l,...,N' (here, N' — N/2 for even N, while 
N' = (N — l)/2 for odd N). Hence, the normal modes of 
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the chain are symmetric (even) or antisymmetric (odd) 
under reflection with respect to the center [53, , such 
that 



(n) , (n) 
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qf 1 ^ , y ■ ' , z- , and n labels the mode. Some 



with w 

general properties can be inferred from this simple con- 
sideration. For instance, the even modes of the axial mo- 
tion are characterized by constant length L of the chain, 
since qff — qL^i ■ For the odd modes, on the other hand, 
the center of mass of the chain, which coincides with the 
chain center, does not move. Clearly, the center-of-mass 
mode, which we denote by , is an even mode char- 
acterized by equal displacements at the positions xf^ of 
the chain. This property and the orthogonality between 
the normal modes lead to the relation 
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for all normal modes with n > 1. 

It is remarkable that also the lowest axial odd mode 
(stretch mode) and its frequency can be exactly deter- 
mined. In fact, taking cx x\ and substituting into 



© one finds 
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where we have used relation @. Therefore, the frequency 

of the axial stretch mode q^ is \J%v and its value is inde- 
pendent of the number of ions iV of the chain. This prop- 
erty has been first demonstrated in [2jJ . It has also been 
observed by numerical evaluation of the normal modes 
of chains up to 10 ions 0,|2(|. Analogously, the trans- 
verse stretch mode, which is the highest odd transverse 
excitation, satisfies y^ ,zf^ cx xf with eigenfrequency 
\J v 2 — v 2 , which is also independent of N. 

We remark that the invariance under reflection im- 
poses different boundary conditions than the ones that 
are usually chosen for a crystal with uniform ion dis- 
tribution. In a crystal that is translationally invariant 
even and odd modes are degenerate and one may choose 
periodic boundary conditions [30j . In presence of an ex- 
ternal potential with central symmetry this invariance is 
broken, apart for the mirror symmetry with respect to 
the center. Hence, at the edges the eigenmodes fulfill the 
relation w-jy = ±^JV', where the sign is determined by 
the parity. 



III. THE SECULAR PROBLEM 

The systematic derivation of an analytic solution of 
Eqs. (10181 is a challenging problem, since it requires to 



take systematically into account the position dependent 
coupling constant and the long range interaction. Nev- 
ertheless, in the limit of large number of ions JV > 1 we 
can make some simplifying assumptions. In this limit, 
in fact, the interparticle spacing aL^Xi) = x 



(0) 
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is a smooth function of the position, and it is inversely 
proportional to the density of ions per unit length 20] , 



n L (xi) = l/a L (xi). 



(11) 



The density of charges for unit length can be evaluated by 
applying the Gauss theorem to a continuous distribution 
of charges, that is assumed to be uniformly distributed 
in an elongated ellypsoid. The resulting one dimensional 
density is [24( 



3 N ( x 2 
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which is defined for \x\ < L, where 2L is the length of 
the crystal at equilibrium. The density (|12l) gives a good 
estimate of the charge distribution in the center of the 
chain for N sufficiently large |20j . The length L is eval- 
uated by minimizing the energy of the crystal, and at 
leading order in N fulfills the relation [2£| 



L(Nf 



NloeN. 



(13) 



In the following, we use these quantities to derive an 
approximate solution for the long- and short-wavelength 
modes in the limit of N 3> 1 ions. Furthermore, we 
compare the results of the derivation with the numeri- 
cal calculations, obtained by solving the eigenvalue equa- 
tions KilSIl for a finite number of ions. 



A. The eigenmodes in the long wavelength limit 

We use the ansatz qi(x,t) = e lut qi{x) in Eqs. © 
and define the rescaled positions ^ = x\°'/L and 
the rescaled interparticle distances a(£j) = a^{xi)jL. 
With these definitions, denoting for simplicity of nota- 
tion qi —> qi, Eqs. © take the form 



2 v 2 



k = - 2 /CoE„ 



(0) 



^ (0) l 3 



(Qi-Qj) (14) 



where we have introduced the dimensionless constant 
2Q 2 2 



mu 2 L(N) 3 3iVlogiV' 



(15) 



If the number of ions is large (N 1), for the long- 
wavelength modes one can approximate the chain by a 
continuous distribution of charges. In this limit, £ is a 
continuous variable varying in the interval (—1,1), while 
the displacement qi — q(£,i) is a continuous function, here 
denoted by g(£). Then, Eqs. l(Ti|) take the form 



(co 2 -u 2 )q(0 = 7 ^IC NI[tq(0} 



(16) 
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where 



where 



^ (^^3(9(0-9(0) (17) 



3(9(0-9(O) 



while n (£ ) = 1 — £ 2 is the density of charges normalized 
to 4/3. Equations ((TTjjl and (|TT|) are valid away from the 
edges of the chain and for long-wavelength excitations, 
where the continuum approximation is reasonable. The 
continuum limit for Eqs. Q gives 



(18) 



A similar type of equations is obtained for Zj. It is re- 
markable that the axial trap frequency enters only as 
a prefactor on the right hand side of Eqs. (|16f) . Con- 
sequently the axial cigcnfrcquencies are proportional to 
v. The transverse eigenfrequencies, instead, do not show 
this behaviour, as one can see from Eq. (|18|) : Here, it is 
the quantity lu 2 — v 2 which is proportional to v 2 . The 
proportionality constants depend on the modes and will 
be calculated in what follows within some approxima- 
tions. The results are independent of the charge Q and 
of the mass m. 

According to Eqs. (|16I18|) . the secular problem consists 
of solving the eigenvalue equation 



Jfc,™ (n) (0] = W n) (0 



(19) 



where u/ n )(£) can be the axial or the transverse modes, 
that satisfy the orthogonality relation 



d£n(£)</j (n) (£)*uA m - l (0 = 5. 



(m). 



(20) 



By partial integration Eq. (|17fl can be written as the sum 
of two terms 



I = h 



A/, 



where Iq contains the contributions of the ions around the 
point £, while the term AI is determined by the value of 
the density and of the eigenmode function and their 
derivatives at the end points of the chain. In appendix A 
we derive the explicit form of the two terms and discuss 
their order of magnitude. For the long wavelength modes 
and at the points £ sufficiently far away from the chain 
end-points £ = ±1, we find / « Iq where 



JoE, ««(£)] = [ logo(0 - - ) (n(£K(£) + 2n'(0«/(O) 



+o(o(0), 



(21) 



and a(£) is an infinitesimal quantity of the order 1/N. 
Using fTTft and (fT2*)l in (f2~T|) we obtain 

Io[£Mt)] = - (log N + log ~ + log(l - e) + |) JHO] 



JH0\ = (i - - 4^'(?) 



(23) 



In the limits of validity of Eq. I|21|l and for N sufficiently 
large, such that log TV 3> 1, Eq. i|22|) can be approximated 
by the leading order in log N 



IoiZMO)-- log NJ[w(0] 
and the eigenvalue equation (|19f) reduces to 



J>W(£)] = A, 



(0 



(24) 



(25) 



where A„ = — log NX n . Equation 125|) is the differential 
equation fulfilled by the Jacobi polynomials P l ' (x) at 
the eigenvalues [3lJ 



A„ = -£(£ + 3) 



(26) 



with £ = 0,1,... and n = £+1. After substitution of (231- 
I26|l into i|19|) , the eigenfrequencies of the axial excitations 
are found from (|15I16|) and take the form 



Jac 



n(n + 1) 



(27) 



Analogously, the eigenfrequencies of the transverse 
modes are obtained from Eq. Ijl8(l with (flfjf) . resulting 



_L Jac 



(n- l)(» + 2) 



(28) 



It should be noted that the solutions of (|23|l for £ = 
and £ = 1 are exact solutions of the original problem JT2J: 
The corresponding eigenmodes w^^x) =constant (cen- 
ter of mass) and w^ 2 \x) oc x (stretch mode), which are 
the continuum limits of the eigenmodes we have found 
for the discrete case, are in fact Jacobi polynomials. The 
result for the center of mass is obvious. The exact re- 
sult for the stretch mode can be understood, noting that 
P 1 1 ' 1 (0 has only one node, whose position coincides with 
the center of the chain £ = and thus with the symmetry 
center for reflections. Hence, its position is independent 
of the number of ions in the chain, and in particular it is 
independent of whether the ions distribution is discrete 
or continuous. 

It is remarkable that the dispersion relation l|27|l co- 
incides with a specific one dimensional limit of a three- 
dimensional continuum mean-field theory, like the one de- 
veloped in |32, although there is no obvious justification 
for this. The two limiting cases, the uniform spheroidal 
fluid of 32] and the case of N strongly-coupled oscillators 
investigated in this work, seem to provide the same axial 
eigenfrequencies in the long wavelength regime and in the 
limit N 3> 1. This result is intriguing, especially if put 
in connection with theory of cold gases in low dimen- 
sions, where different dispersion relations are obtained 
depending on the assumption on the type of mean-field 
interaction |33t l37| . and will be object of future investi- 



(22)gations. 
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FIG. 1: (a) Transverse and (b) axial spectrum of eigenfre- 
quencies (in units of v) for a chain of N = 1000 ions and 
with v t — 400ia The solid line corresponds to the numeri- 
cal solution of Eqs. I6I8H with ||3J. The dashed line shows (a) 
oj„ Jac , (b) wJl Jac : These curves have been truncated, as they 
do not correctly reproduce the short wavelength eigenmodes. 
The dash-dotted line gives the spectra evaluated using the 
method of Dyson |3^] implemented in subsection IHIBI 



Figure ^ presents the comparison between the spectrum 
of eigenfrequencies, obtained by numerically diagonal- 
izing the matrix corresponding to Eqs. iKiiSI) for 1000 
ions, and the formula 1)27(1 and l(28|l derived for the axial 
and the transverse modes. Figure GJa) exhibits the part 
of the spectrum with the long- wavelength axial modes: 
Here, one sees that Eq. 1(2 7JI approximates well the lowest 
part of the axial spectrum, where the limit of continuous 
charges distribution is reasonable. The short-wavelength 
modes are not reproduced correctly by (J2ZJ- These are 
better evaluated by using a more proper approximation 
for this regime, and will be discussed in the next section. 

We remark that, apart for the first two eigenmodes, the 
Jacobi polynomials describe the eigenmode excitation at 
leading order in log N and near the center of the chain, 
where the interparticle separation is of order 1 /N and the 
distribution of charges can be treated as a continuum for 
sufficiently long wavelengths. The continuum approxi- 
mation fails at the edges, where the interparticle spacing 
is significantly larger and Eq. I|12l) is not meaningful. In 
particular, Eq. lfH?)l gives the upper bound for the chain 
length, which would be obtained in the limit of N ^> 1 
particles. Hence, a reasonable boundary condition is to 
assume that the eigemodes and their derivatives vanish 
at x = ±L, where there are no charges and hence the 




FIG. 2: (a) Long-wavelength excitations and (b) Short- 
wavelength excitations of the axial spectrum of eigenfrequen- 
cies. Same notation and parameters as in Fig. 




FIG. 3: (a) 8uj„/un and (b) Suj„/Auj„ as a function of n, 
where Suo n = cjJl' Iac — and AwJl = — The fre- 

quencies ujn are obtained by solving numerically @ with Q • 
From top to bottom: N = 10, 50, 200, 500, 1000. 



energy density is zero. The solution ((23(1 taken at the 
center of the chain neglects the charges at the edges, on 
the basis of the observation that there the number of ions 
is much smaller than at the center, and their contribution 
to the integral l|17|l can therefore be neglected. 

The evaluation of the correction to the results 1(27(1 
and 1(28(1 should be done in perturbation theory in the 
parameter l/log./V, following an analogous procedure to 
the one applied in [2(| for evaluating the correction to the 
density of ions ((12(1 and to the equilibrium length l(13|) . In 
practice, this expansion has a very slow convergence, and 
does not allow for a simple analytical expression. Never- 
theless, the comparison with the spectrum evaluated nu- 
merically, by solving equations J3J and ((618(1 . shows that 
Eq. 1(271) and ((28(1 give already a good estimate of the 
eigenfrequencies for a chain of 10 ions, as it can be seen 
in Fig.[2Ja), where the relative deviation of the frequency 
given by Eq. ((27(1 from the numerical result is plotted for 
chains with different number of ions. This agreement is 
in general valid, respectively, for the axial low modes and 
the transverse high modes, and exhibits a very slow im- 
provement as the number of ions in the chain increases, 
due to the slow convergence of the 1 / log N expansion. 
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B. The density of states in the short- wavelength 
limit 



Simple physical considerations show that the short- 
wavelength eigenmodes are characterized by relatively 
large displacements of the ions around the center of the 
chain, while the ions at the edges nearly do not move. In 
fact, the interparticle distance is minimal in the middle 
of the chain, and for N ^> 1 it is of order 1/N, while it 
is consistently larger at the edges. Hence, a wave cannot 
propagate in a region where the interparticle spacing is 
larger than the wavelength. Here, we can make some sim- 
plifying assumptions in solving the eigenvalue equations 
(KilSt for the short wavelength modes. In fact, in the cen- 
ter of the chain we expect that the relevant contributions 
to the force on an ion originates from the neighbouring 
charges, as nearby groups of charges move in opposite di- 
rections, resulting in forces that mutually cancel. By this 
hypothesis, in 10181) we keep only the nearest neighbour 
interaction, so that the equations to solve are 



q't + v 2 qi = A[qi,q l±1 } 

Hi + vlyi = --A[vi,yi±i] 



Zi + v t Zi 



--A[z i: Zi±i] 



(29) 
(30) 

(31) 



where the tridiagonal matrix A is defined by its action 
on the vector (wi, . . . , ii>i-i, Wi, u>j+i, . . . , u>n) through 



A[ 



Wi,W i± i\ 



Ai(w 



i) + Aj_i(wi_i - Wi) (32) 



and Aj = Kij + i/m where Kij is given in JSJ. The matrix 
A is symmetric, and we denote its characteristic frequen- 
cies by — uj 2 . Following the derivation of Dyson [24]], that 
is summarized in appendix B, the density of states D{Cj 2 ) 
for N ^S> 1 is found from the characteristic function Q(0) 
according to 



D(l/z) = -z 2 Rc 



1 ,. dCl. 
— hm — [—z 
in e->o dd 



(33) 



while Q(6) is explicitly evaluated by using the properties 
of antisymmetric matrices and takes the form 



2JV-1 

^ = iv E iog[i+c(fi,i)]. 



3=1 



Here, £(#, j) is the infinite continued fraction, 



Wi) 



6Ai 



6Ai 



(34) 



(35) 



i+- 



and -Azj-i = A2i = Aj. For N 1 and around the center 
of the chain we may assume A(£) to be a slowly- varying 
function of the position £, such that A i+1 = Aj + <5Aj and 
SAi/Ai <C 1. This allows to evaluate explicitly C($>i) a t 
first order in SAi. For this purpose, we define Aj_i = 



Aj = A, and A J+ i = Aj+2 = A + SA, where SA is the 
first order variation. Furthermore, we denote £(0, j) — 
C and assume ((0,j + 2) = C + where <5£ is a first 
order variation. We substitute these quantities into l|35|) 
keeping only the terms up to first-order and look for a 
consistent solution. The resulting equation is 



(5C(2C + 1) = OS A 
which is integrated to 

C 2 + C - OA = 0. 



(36) 



(37) 



Here, we have taken the integration constant to be zero 
since at the boundaries of the chain A — > 0. The resulting 
solution has the form 



C(o,0 



1 r 



y/l + 40A(Lf) - 1 



(38) 



leading to the characteristic function 



N-l 



,;=i 



3 / d£ n(0 log 



i(Vl + 40A(iO + l) 



(39) 



where we have used the rescaled variable £ and the fact 
that the integrand is even in the interval (—1, 1). Equa- 
tion l|39l) corresponds to the continuum limit of the dis- 
crete summation in Eq. i|34[) . and it is valid away from 
the edges for N > 1. Substituting (J25J) into |J33) we ob- 
tain the equation for the density of states as a function 
of the physical parameters, 



D(l/z) 



(iz 



v /-l + 4zA(L£) 



where f(z) = y/l - (l/4A ;z) 1/3 , while 
A(LO = Ao(l-C 2 ) 3 , 

with 

/Co 



A ° ^ a« = 0)3 
at leading order in log AT. 



2 9 iV 2 2 

V ~ — ; —V 



32 log AT 



(40) 



(41) 



(42) 



Equation l|40(l with Cj 2 = 1/z gives the density of states 
of the short wavelength modes. The corresponding spec- 
trum is shown in Fig.^ The deviations from the numeri- 
cal result are due to the assumption of nearest-neighbour 
coupling, and are small in the short-wavelength part of 
the spectrum, showing that Eq. (|4t)(l provides a good ap- 
proximation in this regime. In particular, this result al- 
lows to evaluate explicitly the value of the maximal axial 
frequency wJLax, the minimal transverse frequency ^ in , 
and the spectrum of the eigenfrequencies in their neigh- 
bourhood. The maximal axial frequency and the minimal 
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transverse frequency are found from the maximal value 
of z for which the integrand in (|4U|) is real, corresponding 
to f(z) — 0. For larger values of z the density of states 
vanishes. The corresponding eigenvalue of the matrix A 
is Q 2 — 4A . From I|29I31[) one finds wJLx = V" 2 + 4A 
and = \J v 2 — 2A . Therefore, the largest value of 
the axial frequency is determined by the largest value of 
the spring constant, that is the value of the spring con- 
stant at the center of the chain, and at leading order in 
logA^ is 




(43) 



Analogously, the smallest value of the transverse modes 
frequency is 



9 



iV2 



16 logN 



(44) 



From Eq. l|4*4*|l one sees that u>^ in can vanish for certain 
values of is, v t , and N. We denote by 



9 N 2 
16 V log A" 



(45) 



the value of the transverse frequency, such that for v t < 
is t CT the linear chain is unstable with respect to excita- 
tions of the transverse vibrations. Using the notation 
introduced in [22] we define the critical value of the as- 
pect ratio between the trap frequencies a cr — v 2 jv 2 cr . It 
takes the value 



16 log AT 
~9~ N 2 



(46) 



which fixes the condition on v t for the chain stability 
according to the inequality v 2 > a~^v 2 , and it is in 
agreement with the analytical estimate in |24j . which 
was obtained under different requirements. In Fig. 0] 
we compare the result (I46f) with the relation cN^ 1 , with 
/3i = —1.73 and c = 2.53, obtained in 22] by fitting 
points calculated with molecular dynamics simulations, 
and verified experimentally in |2^| . Our result reproduces 
approximately this curve. At a cr the crystal undergoes a 
structural transition from a linear string to a zig-zag con- 
figuration, as it has been observed in [6j. It is interesting 
to ask whether this structural transition can be consid- 
ered as a phase transition, as discussed in [23, [U 13 • 
A systematic investigation in this direction requires a 
proper definition of the thermodynamic limit for this kind 
of system. A natural thermodynamic limit, that will be 
discussed in the next section, is one where as N — > oo, 
the axial frequency v — > so that the density in the cen- 
ter is fixed. From (|12|) and l|13fl this requires that the 
ratio v 2 N 2 /log AT is kept constant. From l|45|l it is clear 
that in this limit the critical transverse frequency tends 
to a well defined value. The exploration of the proper- 
ties of this transition in the thermodynamic limit will be 




1000 



FIG. 4: a CI as a function of the number of ions. The solid 
curve gives the result l|460 . The dashed curve is a fit ac- 
cording to the function cA -1 ' 73 , with c = 2.53, as calculated 
numerically in [22^ (see also |23jp. 



object of further studies. Particularly interesting is the 
comparison with standard phase transitions |35| . 



C. The phonon-like approximation 

It is natural to introduce a phonon-like approximation 
for the eigenmodes of Eqs. KilSI) . In this approximation, 
a phonon-like solution is superimposed by a slowly vary- 
ing amplitude, that takes into account the slow variation 
of the coupling strength Kij of Eq. (J5J as a function 
of both i and j. This approximation, that is outlined 
in appendix C, is reasonable for a relatively wide range 
of long wavelength excitations compared to the Jacobi 
polynomials solution, discussed in Sec. lIII Al fsee Fig- EJ - 
It is inferior compared to the Jacobi polynomials in the 
very long wavelenght limit, and it is a bad approximation 
for the short wavelength regime. Therefore, the phonon- 
like approximation, that is natural in condensed mat- 
ter physics, is not an asymptotic approximation in the 
thermodynamic limit N — » oo, neither for the long wave- 
length nor for the short wavelength part of the spectrum. 



IV. STATISTICAL MECHANICS 

In this section we use the density of states, that was 
evaluated in the preceding section, in order to derive the 
thermodynamic quantities of a linear crystal of N ions. 
The linear chain is assumed to be in the regime of stabil- 
ity and to be in equilibrium with a thermal bath at tem- 
perature T. The oscillations around the classical equilib- 
rium points are quantized using standard procedures [2(| . 
It should be noted that in this limit the quantum statis- 
tics of the atoms are irrelevant: In fact, the single particle 
wavepacket is much smaller than the interparticle dis- 
tance [3(j . The dynamics of the crystal is thus described 
by 3 N harmonic oscillators of frequencies uj\ , ui^ , where 
the frequencies u>^ are doubly degenerate. It is modeled 
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by the Hamiltonian H, obtained from H after quantiz- 
ing the eigenmodes of H hal . in Q. Here, H = H + H n , 
where H is the ground state energy, 



N 



H = V + -^ 



(-1 



2ui 



while H n describes the contribution of the collective ex- 
citations 



N 



N 



n=l 



where Nil, N^ y , N„ z are the operators counting the 
number of excitations. The term Vq corresponds to the 
classical minimum energy of the Coulomb crystal. For an 
infinite chain, it is obtained by minimizing the classical 
energy with respect to the length L using the density of 
charges in Eq. I|12|l . and it is evaluated to be [2(| 



V 



10 



mv 2 NL{Nf 



(47) 



where L(N) is given in 

We remark that we investigate the thermodynamic 
quantities for crystals characterized by a finite number 
of particles N and finite axial frequencies v, i.e. crystals 
of finite size, that may be close to experimental situa- 
tions. It is however instructive to consider the definition 
of the thermodynamic limit for this kind of system char- 
acterized by strong correlations, where the effect of the 
charges at the edges cannot be neglected a priori in evalu- 
ating the statistical properties. Here, the thermodynamic 
limit can be appropriately defined by assuming constant 
interparticle spacing (thus constant linear density) at the 
center of the chain x — 0, namely requiring <2l(0) to be 
constant. Denoting by ao = ai,(0), it scales as 



«o = 9 



V vN 



2/3 



where g = (3Q 2 /m) 1 / 3 . This requirement corresponds to 
a vanishing axial frequency, according to v ~ ^/log N /N, 
as N — > oo. With this definition, in the thermody- 
namic limit the maximum axial frequency (|43|> and the 
critical transverse frequency (|45|l are independent of N 

and of v, taking the values wJLax = 3(g/2a ) 3 / 2 and 
z^.cr = 3/4(g/ao) 3 / 2 . In the following, we derive the ther- 
modynamic quantities for ion chains of finite size char- 
acterized by a fixed and finite value of the number of 
particles N and of the axial frequency v, and discuss 
how these quantities behave when we take the specific 
thermodynamic limit that was defined above. 

Assuming thermal equilibrium with the bath, the state 
of the system at constant number of ions N is described 
by the density matrix of the canonical ensemble 



p = — exp (~(3H) 



(48) 



where /3 = 1/ksT and Z is the partition function, 
Z = exp {-0H O ) Y[[l - exp(-/3fkv r , 



which determines the Hclmoltz free energy 

F = -k B T log Z. 

We identify the thermodynamic variables with T, the 
temperature, N, the number of atoms, and v, the axial 
trap frequency, whose variation corresponds to a vari- 
ation in the length of the chain lj. These are not a 
complete set of thermodynamic variables, but they fully 
determine the state of the crystal for the thermodynamic 
quantities we investigate in the following. In particular, 
we take v t as constant and assume wJLax <§C ^min' i- e - that 
there is a large gap between axial and transverse excita- 
tions. In this limit, we consider temperatures T such 
that UbT <SC ?kt>^ in . In this regime the transverse modes 
can be considered frozen, hence the contribution due to 
their excitations to the crystal's thermodynamic proper- 
ties can be neglected, and the dynamics of the system 
is one-dimensional. The thermal energy of the crystal is 
given by 



C/th = (Hn) = E 



hoj'J, 



n exp(/3hu} n ) - 1 
The heat capacity C a — dUth/dT\„.N is 



(49) 



° a ~ dT ^ 



hu>. 



exp(phujn) — 1 



(50) 



The behaviour at high temperatures, such that fc^T 3> 
Swlax (but fc^Tc hcj^ in ), is given by the Dulong-Petit 
law C a = Nks as is clear from (|50|) . On the other hand, 
for ksT <C hp all modes are frozen and the energy of the 
chain is the zero point energy Ho. For large number of 
particles N 3> 1 we can approximate the sum in l|50(l by 
the integral 



C a 



d_ 

df 



exp(^hu) — 1 



where g(u>) = dn/doj is the density of states. In partic- 
ular, at temperatures such that the contribution of the 
long-wavelength excitations to the sum is predominant 
and is given by l|27ll . the resulting density of states is 



4lu 



low T 



and the heat capacity is given by the integral 



C a \ 



V2k% d 



B W rp2 



low T hv Q T 



dx- 



1 V^ 2 + 4/8 



(51) 
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FIG. 5: Specific heat c a — C a /N, in units of fcs, as a function 
of the temperature T, in units of fiLJm&x/kB. Here, N = 1000 
and wJLax -C u^ in . This would correspond to an experiment 
with v = 1-7TX 1 kHz and v t — 2nx4 MHz. The inset shows the 
low-temperature behaviour, the dashed line is the estimated 
slope, with c = y%C t (/£)k%/%v (see text). 



where we have defined xq = (3Tiv. Hence, the integrand 
and the integration limits depend on the temperature 
through xq. In this regime and for /c^T 3> fiv we can 
set x Q ~ in l|ET|l. and recover the result C a = cT, with 
c = \/%^{2)k^/hv, where £(2) is the Riemann's zeta func- 
tion. Therefore, for the considered regime the heat capac- 
ity is proportional to the temperature, which is a charac- 
teristic behaviour encountered in a one-dimensional De- 
bye crystal [3(j. In Fig. the specific heat c a — C a /N 
for N = 1000 ions is plotted as a function of the temper- 
ature T. In the inset, the low temperatures behaviour 
is shown, and compared with the curve cT/N estimated 
above. The figure shows that the evaluated behaviour, 
valid in the asymptotic limit of an infinite number of ions, 
provides a reasonable description of the specific heat at 
low temperatures. 

It is remarkable that the heat capacity in (|51|l scales 
like C ~ 1/v. The specific heat per particle c a — C a /N 
behaves thus like c a ~ 1/Nv at low temperatures, and in 
the thermodynamic limit it vanishes as c a ~ l/yflogN. 
It thus depends on the number of ions, and this is a mani- 
festation of the deviation from extensivity of the system's 
behaviour. Note that the relative energy fluctuations are 
of the order (\/\og N /TV) 1 / 2 , therefore the usual equiva- 
lence of ensembles holds. 



The pressure P in the axial direction is defined as 



P 



dF 
dL 



T,N 



(52) 



and it is the variation of the free energy with the length of 
the string at constant N and T. Under these conditions, 
the length of the crystal is varied by changing the axial 
trap frequency v, according to 



dL 
dv 



T,N 



2L 
~3v~ 



(53) 



as obtained from ((13(1 . Substituting the explicit form of 
the free energy into 1)520. we find 



where 



P, = 



dH 



dL 



P = P + Pt, 



Vq 3h x ^ 

= — H > 

T.N L 4L ^— ' 



2u± 1 - 



,X2 



(54) 



is the pressure at zero-temperature, and Pt gives the 
contribution of the excitations, 



Pt - ^U th . 



(55) 



In deriving l|54|) and (|55|) we have used dVo/dL = —Vq/L 
that results from 1)47(1 . and the relations 



ov 



d_ 

dv 
d_ 

dv 



V 
V 



n 



1 - 



,,-L 



that are obtained from the functional dependence of ujn 
and ui^ on v, as can be extracted from (|16(l and 118(1 . 
Note that a variation of the axial trap frequency implies 
also a variation of the transverse eigcnfrcqucncics, which 
give a contribution of opposite sign to the pressure, as it 
is obvious from l(54|) . However, the contribution due to 
the quantum mechanical zero-point energy is very small 
compared to the classical term Vq/L. Therefore, Po is 
dominated by Vq/L and, using ((47|) . in the thermody- 
namic limit it scales as Po ~ logA^. The term Pt de- 
pends on the temperature. For low temperatures, such 
that UbT ^> hv, it scales as Pt ~ 1/Lv ~ 1/^log N. 
At high temperatures, in the Dulong-Petit regime, Pt 
does neither depend on N nor on v. In the regime where 
the chain is thermally stable, which we consider here, 
Vq » U t h giving P w P . Thus, the pressure is dom- 
inated by the zero-temperature contribution and in the 
thermodynamic limit P pa Po ~ log N. A useful relation 
for the following discussion is 



dP 
df 



L.N 2L 



(56) 



The isothermal compressibility kt is evaluated from 
the pressure according to 



1 

kt 



B 



dP 
"dL 



T.N 



(57) 



where B is the bulk modulus. Using 1(54(1 . 1(55(1 in 1(57(1 

we find 



kt 



(58) 
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In the thermodynamic limit the bulk modulus B is dom- 
inated by the zero-temperature contribution —LdPo/dL, 
that in turn is dominated by the term —LdVo/dL ~ 
Vq/L ~ log AT. Therefore, B ~ log A" and the compress- 
ibility kt vanishes as 1/ log AT. 

The coefficient of thermal expansion q?t can be evalu- 
ated from the knowledge of kt and C a according to [3(j 

IdL 1 dP/dT\ L 3 

L dT p,n LoP/dL\T 2L 

where we have used (I56|) . Since the compressibility is 
dominated by the zero temperature term, the behaviour 
of the coefficient of thermal expansion q-t as a function 
of T is determined by the heat capacity: Linear depen- 
dence at low temperatures, and saturation at high tem- 
peratures. At low temperatures C a ~ l/i/, and as the 
thermodynamic limit is approached ar ~ (logA^)~ 3 / 2 . 
At higher temperatures, when the heat capacity mani- 
fests the Dulong-Petit behaviour C a = Nks, the coeffi- 
cient of thermal expansion vanishes like cut ~ 1/logiV. 
In Fig. the coefficient of thermal expansion for a chain 
of 1000 ions is presented. Calculations made with dif- 
ferent numbers of ions, taking the trap frequency v such 
that the linear density at the center of the chain remains 
constant, show that clt decreases with N. The numer- 
ical results are consistent with the behaviour we expect 
in different temperature regimes, according to the above 
considerations. In particular, for finite A~ it is signif- 
icantly different from zero. This is in contrast to the 
behaviour found for uniform harmonic solids, where the 
coefficient of thermal expansion vanishes 30] . 

The thermodynamic quantities of a linear string of 
charges confined by a harmonic trap are affected by the 
way the thermodynamic limit is taken. Nevertheless, the 
behaviour of the system is intrinsically non-extensive. 
The non-extensitivity is due to the strong correlation 
and the dimensionality of the crystal, which determine 
a regime where the correlation energy, associated with 
the discreteness of the individual charges, cannot be ne- 
glected |2(j . It manifests in particular in the log N terms 
appearing in the thermodynamic quantities. A represen- 
tative example is the dependence of the specific heat per 
particle on N. 

We finally remark on the thermal stability of the chain. 
The derivation presented in this section, in fact, relies on 
the assumption that the thermal excitations do not affect 
the stability of the system and thus the validity of the 
physical model we are considering, namely of ions oscil- 
lating around their equilibrium positions. This condition 
is equivalent to the statement that the thermal energy at 
the considered temperatures is considerably lower than 
the equilibrium energy, and the displacements are much 
smaller than the respective interparticle distance. Hence, 
this condition is valid when 

O 2 

— log A^ > k B T. 
a 

This relation is amply satisfied for the parameters of the 




FIG. 6: Coefficient of thermal expansion olt, in arbitrary 
units, as a function of T, in units of ftajm ax /fcs. Here, 
N = 1000 and wJLax <C ti),J; in . This would correspond to an 
experiment with v = In x 1 kHz and v t = 2tt x 4 MHz. For 
these parameters and Berillium atoms, T is in units of the 
Debye temperature Od, with Od = ftuiai/fe ~ 20/iK, and 
«t in units of 10 _5 /iK _1 . 

numerical evaluation of the thermodynamic quantities 
considered in this section. 



V. SUMMARY AND CONCLUSIONS 

We have investigated the dynamics of a linear crys- 
tal excitations with Coulomb interaction and in presence 
of an external potential. We have derived an analyti- 
cal formula for the density of states at long- and short- 
wavelengths. In the long- wavelength part of the spec- 
trum, we have calculated analytically the eigenmodes and 
eigenfrequencies. The eigenmodes and eigenfrequencies 
for the center-of-mass and the first excitation of the axial 
and of the transverse motion are exact and independent 
of the number of ions. Apart for these modes, the re- 
sults we derive are valid in the limit of infinite number of 
ions. Nevertheless, they already give a good description 
of the spectrum of excitations of chains of dozens of ions. 
Using our results we study the statistical mechanics and 
thermodynamics of the linear chain. 

Our derivation allows to find an analytical formula for 
the critical transverse frequency required for determin- 
ing the stability of the linear chain. It agrees with the 
analytical estimate by |24|, that was obtained under dif- 
ferent assumptions, and it is consistent with the formula 
fitted from numerical data and verified experimen- 
tally • It was suggested that this instability of the ion 
chain, resulting in a transition from a linear to a zig-zag 
equilibrium configuration, can be treated as a phase tran- 
sition [22H23I]. In future works we will explore, using the 
formulation developed in this work, whether the thermo- 
dynamic quantities exhibit singularities characteristic of 
phase transitions j3^. This system differs from systems 
that are traditionally studied in the framework of statis- 
tical physics, since it is not extensive. In particular the 
specific heat per particle depends on the number of ions. 



11 



The results presented in this work show a statisti- 
cal mechanics approach applied to a strongly-correlated 
mesoscopic system. They contribute to the on-going re- 
search on low-dimensional cold gases and may be 
relevant to studies of the quantum dynamics of many- 
particle Coulomb systems like ion crystals in storage 
rings and cold neutral plasmas [39l Eof , which at suffi- 
ciently low temperatures are predicted to crystallize |4l| . 
The connection with Wigner crystals, where the quantum 
statistics may play a relevant role [42|, will be explored. 

Moreover, the spectra of excitations here evaluated are 
relevant for the implementations of quantum logic with 
ion traps, where the knowledge of the long- wavelength 
modes is important for the realization of logic gates Ho . 
OEM 113 ■ as well as realization of solid-state models [43| . 
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APPENDIX A: SECULAR EQUATION OF THE 
LONG-WAVELENGTH EIGENMODES 

In this appendix we outline the derivation of the dif- 
ferential equation I|22I23[I from Eq. (|17fl . and estimate 
the correction to the solution we find. By integration by 
parts, Eq. I|17|) takes the form 



1 (m-i) /«(!) 



where 

h = 



1 



2a(£) 2 
1 



I = I + L + h 



(/e(£-<0 + /*(£ + <*)) 



(Al) 



2 



2a(0 

^loga(0 (jf + if (£ + «) 



2 V(i + 2 (i-O' 



+\ (if iog(i + + /| 2) (i) iog(i - 0) 



d£'iog(£-Oif(0 



d£'iog(r-o/i 3) (n 



(A3) 



and = n(£')(w(£) - u>(£')). Clearly the integral I 

results of the contribution of the ions around the ion at 
£. For long-wavelength modes and away from the end- 
points £ = ±1, the interparticle spacing a(£) scales as 
1/N and it can be treated as an infinitesimal quantity. 
Expanding IjAlfl in a(£) and keeping the leading order, 
we obtain Eq. I|21|l . The integral (|A2|) contains the con- 
tribution of the end points of the chain. Also, the inte- 
gral (|A3f) is dominated by the end points, as can be seen 
by integrating repeatedly by parts, taking the indefinite 
integral of logx, and making the reasonable assumption 
that f(x) is an analytic function on the interval (—1, 1). 

By imposing the requirement that the density of ions 
and its derivatives vanish at the end-points, while 
and its derivatives are bounded, the term AI = I\ + 12 
is much smaller than Iq. In particular, the contributions 
at the end-points vanish, whereas the contributions to I2 
from the ions around the point £ are of order a(£) ~ 1/7V. 



APPENDIX B: CHARACTERISTIC FUNCTION 
FOR THE SHORT-WAVELENGTH SPECTRUM 

In this appendix we discuss the evaluation of the char- 
acteristic function for the eigenvalue problem in Eq. I|29() 
following the treatment developed by Dyson |34|. The 
procedure can be immediately generalized to the eigen- 
values of the transverse motion. Substituting qj(t) — 
J duie lut qj(uj)/2ir into l|29|l and replacing qi(uj) with g^, 
we obtain 

- {uj 2 - v 2 )q 3 = Aj(q j+ i - qj) + Aj^qj-x - q 3 ) (Bl) 

The roots — tl> 2 of the matrix A are related to the eigen- 
frequencies ui by 



(A2) 



u> = \/ uj 1 + v 2 

These roots can be found by solving the second-order 
differential equation 

% = A j fe+i + A i-i fe-i - * ) ( B2 ) 

where qj(t) = J dQe iGjt q J (Q)/2ir. Following 0, this 
problem reduces to the diagonalization of an antisym- 
metric matrix. In what follows we review the funda- 
mental steps. We define the variables si, . . . , sn-i, such 



that Si 



1i, 



After substituting into 
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Eq. I|B2|) and integrating, we get the first-order differ- 
ential equation q~j = 



duce the variables u\ 
u 2] = s 



— y/ kj-\Sj-\. We intro- 
,U2N-i, such that it2j_i = <jj, 



and the new matrix elements Aj, such that 



A 2j - A 2i _i 



A,-. With these definition, we have now a 



set of 2N — 1 first-order differential equations 



which are now defined by the (2N — 1) x (2N — 1) 
antisymmetric matrix A, whose elements arc 

Aj+ij = — Ajj+i = \\jhj. The characteristic fre- 
quencies ujj of the chain are the square roots of the 
characteristic roots of —A, of which one is and the 
other 2 N — 2 come in pairs ±£jj . 

Using the matrix A, in the limit JV>1 the spectrum 
of characteristic frequencies M(lo 2 ), that is defined as the 
proportion of roots with cD| < a) 2 , and the corresponding 
density D(u> 2 ) = dM(fi)/dfi\ li= Qa, can be derived from 
the characteristic function £1(9), which is defined as |34j 



2N-1 



0(0) 



lim 

JV^oo 2N 



l — E iog(i 



OCo 2 ) 



(B3) 



where 9 is a complex variable. The density of character- 
istic frequencies -D(/i) and the spectrum M(fi) are found 
by using the properties of the analytic continuation of 
fl(0), 



Re 



— lim il(~z + ie) 



djuD(/x) = 1 - M(l/z) 



from which Eq. I|33|) is obtained. 



APPENDIX C: AN ANSATZ FOR THE SECULAR 
EQUATION: SLOWLY VARYING ENVELOPE 

In this appendix we discuss an ansatz based on the 
phononic solution for a translationally invariant crystal, 
where we assume that the envelope, superimposed on the 
amplitude and the phase of the phononic solution, varies 
slowly as a function of the position. This ansatz is sup- 
posed to be valid for long-wavelength excitations and 
chains with N 3> 1 atoms. 

The exact solution of in the limit of uniform 

charge distribution with constant spacing a, Kij/m = 



j\ A with k = 2Q 2 /a A , is |3C 



1j 



^i(kja-ujt) 



(Cf) 



Here, to is the eigenfrequency, k the wavevector, and the 
boundary conditions for the uniform chain are fixed so 
that qN = qo, giving e lkNa = 1. Hence, ka = 2ttI/N with 
I = 0, 1, . . . , N — 1. The dispersion relation is obtained 



by substituting (|C1(I into © with Kij/m = Kj 



-J 



If 



the interaction is determined by the nearest-neighbours, 
while the interaction with the other ions is neglected, 
then the eigenfrequencies are given by 



= 2V^ 



ka 



sm ■ 



(C2) 



Note that u> depends on ka, which takes the values 
ka = 2ttI/N. 

We now construct from this solution an ansatz for 
the inhomogeneous chain, in the limit of slowly- varying 
spring constant Kjj+i over the wavelength of the propa- 
gating perturbation. We assume thus the local solution: 



1j = A i e 



i(kja-u)t) 



(C3) 



where ka is a constant, i.e. we have assumed kjXj = 
jka. The non-linear variation of the phase with the site 
is included in Aj. The ansatz l|C3|l assumes that it is 
possible to write the solution as the product of a slowly 
varying amplitude and phase, represented by the factor 
Aj, and a fast oscillating part, that is the exponential 
term in (|C1I) . The validity of this assumption is checked 
later on. 

By substituting (|C3I) into 10, keeping only the nearest 
neighbour interaction, we obtain: 



(co 2 v 2 )A d = Kj (A, - A J+1 e ika ) 



(C4) 



— ika 



ie 



where Kj = Kjj+i/m. We make use of the fact that Kj 
varies slowly with the position, which allows to expand 
Aj and Kj in Sj , of order unity, according to 



A 



Kj 
Kj-l 



A±5A 

k + Ak + 5k 3 /2 
k + Ak — SKj/2 



(C5) 
(C6) 
(C7) 



where 8 A — |jy<$7 with Sj = ±1. The spring constant 



has been divided into three contributions: k is constant 
and fulfills l|C2(l . Ak is also constant and 6k j — K'Sj, 

Ak ^> Skh. The chosen 



with 



^r 2 -, such that k 



variation of Kj reflects the symmetry under reflection of 
the crystal, such that the ion at the center is character- 
ized by equal couplings on both sides (which implies the 
condition 5k — 0). The term Ak does not depend on j, 
and it is the correction to k, the spring constant fulfilling 
the dispersion relation for the case of a uniform chain: 



A, 



k. Substituting l|C5IC7|) into (|CH) yields 



AAkA sin 2 



ka 



2ik5A sin ka — 2iA5k sin ka = 



where we have applied the dispersion relation (|C2fl . Since 
there is only one zero-order term in j, it is evident that 
Ak = 0: Thus, there is no zero-order correction to the 
spring constant k of the uniform chain, and its value 
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in this order is accounted for by l|C2|l . The first or- 
der expansion gives thus a differential equation (valid for 
sin ka =/= 0) 



5A _ 5k 
which admits the solution 



The eigenmode has then the form 

1j - 



K \kja-iujt 
e J 



(C8) 



(C9) 



(CIO) 



and the wave vector takes the same values as for a uni- 
form crystal with half-periodic boundary conditions, tak- 
ing into account the symmetry under reflection, 



kNa = nl 



(Cll) 



where k = k(uj) is given by Eq (|C2|) . and 
I = 0, 1, . . . , N — 1. Therefore, this ansatz gives simply 



a variation of the amplitude of the displacement of the 
ion as a function of the local spring constant, but no 
change in the phase, which is the same as the one of a 
uniform crystal. The amplitude is smaller at the center of 
the chain, where the ions are closer and the coupling con- 
stant is larger. Taking kq — Aq as given by (|42|l . i.e. the 
maximum value of the spring constant in the chain, and 
using l|C2|l . we obtain the maximal frequency as given 
in H43|) . A similar argument can be developed for the 
transverse frequency. 

Figure[7|compares the prediction of the slowly-varying 
ansatz with the numerical results and the Jacobi poly- 
nomials results given by Eq. H27|) . Here, it is obvious 
that the Jacobi polynomials provide a better approxima- 
tion for the spectrum at the lowest eigenfrequencies. The 
curve evaluated from Eq. !jC2(l using Eqs. (|C3fl and IIC11I) 
lies close to the spectrum evaluated numerically for a 
larger range of modes, but it does not reproduce its 
asymptotic behaviour for N — > oo. 
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